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Formal series and numerical integrators: 
some history and some new techniques 

Jesus Maria Sanz-Serna and Ander Murua 


Abstract. This paper provides a brief history of B-series and the associated Butcher 
group and presents the new theory of word series and extended word series. B-series 
(Hairer and Wanner 1976) are formal series of functions parameterized by rooted trees. 
They greatly simplify the study of Runge-Kutta schemes and other numerical integrators. 
We examine the problems that led to the introduction of B-series and survey a number of 
more recent developments, including applications outside numerical mathematics. Word 
series (series of functions parameterized by words from an alphabet) provide in some 
cases a very convenient alternative to B-series. Associated with word series is a group Q 
of coefficients with a composition rule simpler than the corresponding rule in the Butcher 
group. From a more mathematical point of view, integrators, like Runge-Kutta schemes, 
that are affine equivariant are represented by elements of the Butcher group, integrators 
that are equivariant with respect to arbitrary changes of variables are represented by 
elements of the word group Q.. 
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1. Introduction 

This contribution presents a survey of the use of formal series in the analysis of 
numerical integrators. The first part of the paper (Section 2) reviews a number 
of developments (the study of the order of consistency of Runge-Kutta methods, 
Butcher’s notions of composition of Runge-Kutta schemes and effective order, etc.) 
that underlie the introduction by Hairer and Wanner [55] of the concept of B-series 
in 1976. B-series are formal series parameterized by rooted trees. We also show 
how B-series and the associated Butcher group have gained prominence in view 
not only of their relevance to the analysis of geometric integrators but also of their 
applications to several mathematical theories not directly related to numerical 
analysis. In the second part of the paper (Sections 3-5) we summarize the recent 
theory of word series and extended word series developed in [57]. These series are 
parameterized by words from an alphabet, rather than by rooted trees, and may 
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be used advantageously to analyze some numerical integrators, including splitting 
algorithms. Word series appeared implicitly [17] and explicitly in [18] . [19] . Ex¬ 
tended word series were introduced in m- While series of differential operators 
parameterized by words (Chen-Fliess series) are very common in control theory 
and dynamical systems (see e.g. [30], [21]) and have also been used in numerical 
analysis (see, among others, [33], [35], [2D]), word series, as defined here, are, like 
B-series, series of functions. 

Word series and extended word series are, when applicable, more convenient 
than B-series. Associated with word series is a group Q of coefficients with a 
composition rule simpler than the corresponding rule in the Butcher group. From 
a more mathematical point of view, integrators, like Runge-Kutta schemes, that 
are affine equivariant are represented by elements of the Butcher group, integrators 
that are equivariant with respect to arbitrary changes of variables are represented 
by elements of the word group Q. 

It should also be mentioned that B-series and word series have been recently 
applied to reduce efficiently dynamical systems to normal forms, to perform high- 
order averaging and to find formal conserved quantities PS], HZ], nsg, HU, EZj. 
The approach via B- and word series also makes it possible to simplify the deriva¬ 
tion of the corresponding error estimates. These developments provide applications 
outside numerical mathematics of tools originally devised to analyze numerical 
methods. 


2. Series based on rooted trees 

This section contains an overview of the use B-series in the analysis of numerical 
integrators. 

2.1. Finding the order of a Runge-Kutta method. A Runge- 
Kutta (RK) method with s stages is specified by s + s 1 2 real coefficients , 

i,j = 1,..., s. When the method is applied with stepsize h to the initial value 
problen]]] 

4 x = f{x), x(0) = x 0 € R d , (1) 

at 

and the approximation x n to the exact solution value x(nh), n = 0,1,..., has been 
computed, the formulas to find the next approximation x n +i are 


A n ,i — %n T l 1 ^ ajj f (X-n.j ), 1 i ^ s, 


3 =1 


( 2 ) 


X n+ l = x n + h ^2 bif(X n ,i ) 


1 Our attention is restricted to deterministic differential equations in Euclidean spaces. We 

do not consider the extension to stochastic differential equations, differential equations on Lie 
groups, differential-algebraic equations, etc. 
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Table 1. Rooted trees with < 4 vertices. For each tree the root is the bottom vertex 


M 

1233 4 4 4 4 

a(u) 

1112 1 2 1 6 

u\ 

1 2 6 3 24 12 8 4 


. i i v : 

: Y 


(the X n> i are the auxiliary internal stages). If / is well behaved and \h\ is small, 
the relations m define x n +i as a function of x n , i.e. x n +i = iph{x n ). The /- 
dependent mapping iph is an approximation to the solution flow (f>h of (JT]) for 
which x((n + 1 )h) = <j>h(x(nh )); for smooth / and an RK method of order is, 
iph(x) — (j>h(x) = 0(h u+1 ) as h —> 0. To determine the value of is for a given 
method it is therefore necessary to write down the expansion of iph{x) in powers of 
h, a job trickier than one may think: it took mathematicians more than fifty years 
[B| to come up with an easy, systematic way of carrying it out. After the work of 
Butcher [3j —following on earlier developments by Gill and Merson— the Taylor 
series for is written, with the help of rooted trees, in the form 


iphix) =x + ^2h n 

n=l ueTn 


a(u) 


,X u (x). 


( 3 ) 


Here: 

• T n denotes the set of all rooted trees with n vertices (see Table [T]) and, for 
each rooted tree u, a(u ) is the cardinal of the group of symmetries of u. 

• F u (the elementary differential associated with u) is an K^-valued function 
of x that depends on © but does not change with the RK coefficients bi, 
ciij. Using the rooted tree as an example, 


T 


j(x) =d xx f(x) f(x),d x f{x)[f{x)] 


where d xx f(x)[-, ■] and d x f(x)[-] respectively denote the second and first 
Frechet derivatives of / evaluated at x. The key point is to observe how 
the structure of the elementary differential mimics that of the corresponding 
rooted tree. 


• c u (the elementary weight associated with u) is a a real number that changes 
with the coefficients bi, aij and is independent of the system © being inte- 
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grated. Taking again V as an example, 


s s s s 

c a = bi ^ a,ij aik ( dki )^; 

i= 1 j = l fc=l ^=1 

the structure of the rooted tree is reflected in the summations. 
For the solution flow (f>h of CD, the expansion is 


(t> h {x) =x+^2h n y 

n =1 u&Tn 


1 1 

cr(u) u\ 


u(*e)' 


(4) 


where u\ (the density of u) is an integer that is easily computed recursively. By 
comparing m and £3} we conclude that an RK method has order > v if and only 
if c u = 1/zi! for each u &T ni n < v. As an illustration, we note that for order > 3 
the coefficients have to satisfy the following set of order conditions 


y ' hi — i, y' bidij 

i ij 


ijk 


bidijdjk g, 


y ' bidijdik ^ . 
ijk 


These equations may be shown to be mutually independent [7]. 


2.2. Composing Runge-Kutta methods. Processing. Butcher [5] 
developed an algebraic theory of RK methods. If Y 1 and represent two RK 
schemes with si and S 2 stages respectively, the composition of mappings oi/j^ 
corresponds to a single step of size 2h of a third RK scheme (with S 1 +S 2 stages). As 
distinct from the /-dependent mapping o , this new scheme is completely 
independent of the system © under consideration and, accordingly, it is possible 
to define an operation of composition between the RK schemes themselves. The 
study of this binary operation is facilitated by considering classes of equivalent RK 
schemes rather than the schemes themselves, where a class consists of a scheme 
and all others that generate the same 1 ^ (note that, for instance, reordering the 
stages of the method changes the value of the atj and but does not lead to an 
essentially different computation, see ED- Furthermore, if y = iph{x) is an RK 
map, then x may be recovered from y by taking a step of size —h from y with 
another RK scheme, the so-called adjoint of the original scheme |43| . This should 
be compared with the situation for the flow, where y = <ph{x) implies x = <j)-h(x). 

In order to see that these developments are of relevance to practical computa¬ 
tion, we discuss briefly the idea of processing also introduced by Butcher [3]. If 
Xh is a near-identity mapping in R 13 and i[>h represents any one-step integrator for 
©, the mapping 

= Xh 1 °iph°Xh 

defines a processed numerical integrator (this is easily interpreted in terms of a 
change of variables, see [32]). For m > 1, 

i’h = °^ho xh) m = Xh 1 °^h 0 Xh ; 
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therefore to advance m steps with the method ifh one preprocesses the initial con¬ 
dition to find Xh(xo), advances m steps with the original method and then post¬ 
process the numerical solution by applying yT 1 . Postprocessing is only performed 
when output is desired, not at every time step. If both %f>h and \h correspond to 
RK methods, then so does fh as discussed above. The idea of processing is useful 
in different scenarios. If \h may be chosen in such a way that iph is more accurate 
than the original if>h, one obtains extra accuracy at the (hopefully small) price of 
having to carry out the processing; iph is said to have effective order v [ 8 ] when 
'iph has order 9 larger than the order v of ifh- As a second possibility, the pro¬ 
cessed integrator may possess, when iph does not, some of the valuable geometric 
properties presented below. 


2.3. B -series. Butcher series (B-series for short) were introduced in [28] as a 
means to systematize the derivation and use of expansions like ([3]) or (J3I) - It is 
convenient to introduce an empty rooted tree 0 with cr( 0 ) = 1 , 0 ! = 1, (x) = x 

and, for each RK method, c@ = 1. If T represents the set of all rooted trees 
(including 0 ), then © and © become respectively: 


iph(x) = V 9 |u| —!!— c u T u (x), 
£ r a{u) 


4>h(x) = Y hH 

u&T 


1 1 

cr(u) u ! 


Bu, (:r) 


(|w| is the number of vertices of u; |0| = 0). The right-hand sides of these equalities 
provide examples of B-series: if S £ R 7 " (i.e. <5 is a mapping that associates with 
each u G T a real number S u ), the corresponding B-series is, by definition, the 
formal serie^l 

B s (x) = V h} u \-^-8 u F u {x). (5) 

Z£r a{u) 

Note that B-series are relative to the system © being studied because the elemen¬ 
tary differentials change with /. 

Obviously the set of all B-series is a vector space. A more important algebraic 
feature is that, if 7 € MT is a family of coefficients with 70 = 1 and S € R 7 ”, then 
the composition B^[B^(xyj is again a B-series B^(x); furthermore, the elements 
( u of the family ( are functions of the elements of <5 and 7 and do not vary with h 
or /. For instance for the rooted tree 


= <*0 7^ + 8 . 7. 7; + 8 : 7j + <5, 7. 7. + <*j 7. + 8 V 7. + ^ 70- (6) 

The rooted tree in the left-hand side has been ‘pruned’ in all possible ways (in¬ 
cluding complete uprooting —> 0 and no pruning —► '})■ in the right-hand side 
8 is evaluated at the rooted tree that remains after pruning and 7 is evaluated at 
the pieces that have been removed. 

As a first example of the use of B-series, we outline the derivation of the RK 
expansion ©• We begin by assuming that the RK solution x n+ i and each in¬ 
ternal stage X n< i are expressed as B-series (evaluated at x n ) with undetermined 

2 Rather than using the normalizing factor 1 /cr(u), the original paper m uses an alternative 
factor. The present normalization simplifies many formulas. 
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coefficients c £ R 7 ", C l £ R 7 ". Those series are then substituted in m and the 
composition rule we have just described is used to find the B-series for 
this procedure yields equations that determine the elementary weights c u . The 
expansion (j4j) may be derived analogously. 

2.4. The Butcher group. If 70 = 1 and Bs[B 7 (x)) = B^(x) as above, we 
may write £ = < 5*7 (recall ( does not change with / or h). The set of all 7 £ R 7 " with 
70 = 1 is a group Qb for the operation *: the so-called Butcher group. The (Taylor 
expansion of the) true solution flow </>/ 1 and (the Taylor expansion of) the RK map 
%/jh correspond to the elements of Qb {1/u!} and {c u } respectively; equivalent RK 
schemes give rise to the same element of Qb- For points in Qb associated with RK 
schemes the operation * obviously corresponds to the composition of RK schemes 
defined by Butcher. 

Assume that the system (JT]) undergoes an affine change of variables x = Mx + c 
and becomes ( d/dt)x = f{x), with f(x) = M~ 1 f(Mx + c). If 7 £ Qb, B 1 {x) = 
MB^tjc) + c, where the notation B indicates that the elementary differentials in 
the series are based on /; in this way B 7 is affine equivariant. In general!! I ? 7 is 
not equivariant with respect to more general changes x = x{%)• 

Taylor series methods, multiderivative RK methods and other one-step inte¬ 
grators x n+ \ = i/Jh(xn) not of the form ([ 2 ]) are also represented by elements of Qb- 
The characterization of the family of integrators that correspond to elements of Qb 
has been obtained only recently [33); note that those integrators necessarily must 
be affine equivariant. 

A forest is a formal juxtaposition of nonempty rooted trees, such as 
or {\?, 1 ,:}; it is understood that here the order is irrelevant (e.g. {.,, 1 ,and 
are the same forest). Also considered is an empty forest that contains 
no rooted trees. Thus the set of forests T is the commutative monoid generated 
by the set of nonempty rooted trees; by taking all formal linear combinations with 
real coefficients of forests we obtain a commutative, associative algebra A. The 
decomposition of each rooted tree that is necessary to compute Bs{B 1 {x)) (for 
instance from d 6 j) ^ is decomposed as 


0<S>\ J + .(g).: + :(g>: + :®.. + i<8>. + V<g). + \ J ®0) 

defines a coproduct that turns A into a Hopf algebra: the Connes-Kreimer algebra 
that appears in renormalization in quantum field theory and elsewhere. Since T 
is a basis for the vector space A, each linear form on A may be identified with 
a mapping £ R^. Each element 7 £ Qb C R 7 " may be extended to a mapping 
£ R^ by defining 7 at a forest as JI u 7 “’ w here the product is extended to all the 
trees in the forest. In this way the elements of Qb may be seen as the characters 
(multiplicative morphisms) of the Hopf algebra A. The operation * defined above 
is the restriction to Qb of the convolution product in the dual of A. In this way 
Butcher’s work anticipated many results that were to be later used in different 
subfields of mathematics [5]. 

3 Of course the B-series for cf>h is equivariant with respect to all changes of variables. 
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2.5. Modified equations. Since the properties of maps (discrete dynamical 
systems) are often more difficult to investigate than those of differential equations 
(continuous dynamical systems), it may make sense, given a one-step integrator 
iph, to seek a differential system ( d/dt)x = fh{x) whose h flow (/h coincides with 
the map i/>h- While it is well known [38] that it is not possible in general to find 
such /ft, for typical integrators and smooth /, one may construct a formal series 

Jh{x) = f i0 \x) + hp\x) + h 2 P\x) + ■■■, 

whose formal h -flow exactly matches the Taylor expansion of /^(x). For instance, 
for Euler’s rule with i/jh{ x ) = x + hf(x ), /ft(x) is found to be 

f(x) -t \dxf(x)[f{x)) + y <9 x /(x) \d x f(x) [/(x)]] + y d xx f{x ) [f(x) ,/(x)] + ■■■ . (7) 

The fact that fh(x) — /(x) = 0(h) indicates that the integrator has order 1; for 
an integrator of order z/, fh(x ) — f(x) = 0(h u ). When the terms of order h^, 
/i = 1 , 2 ,..., and higher are suppressed from 0, one obtains a vector field fj^ ( x ) 
whose h -flow differs from iph in 0{h^ +1 ). Thus Euler’s rule applied to 0 is an 
approximation of order /i to the modified differential system (d/dt)x = /J^(x); 
for /i large and \h\ small the modified system may be expected to provide a very 
accurate description of the behaviour of the numerical solution. The idea of using 
modified systems is very old, see e.g. the references in [24], and, as we shall see 
later, has gained prominence with the growing interest in geometric integration. 

Note that, for /ft(x) in 0, hfh(x) is a B-series. In fact [25], |TD]., each B 1 {x) 
with 7 & Qb coincides with the h -flow of a uniquely defined vector field /ft(x) 
where hfh(x) is a B-series Bp{x) with /3g = 0. If we denote by gs the set of such 
/3’s, the mapping 7 1 —> ft is a bijection from Qb onto qb ■ Thus B-series integrators 
may be handled either by using the coefficients { 7 ^} of the B-series for ?/>ft(x) or 
by the coefficients {/ 3 U } for the corresponding hfh(x). The second option is often 
advantageous because, while g is a linear space, Qb is not; this implies that in 
many situations properties that are nonlinear when expressed in terms of the "f u 
become linear for the corresponding /?„. Formally Qb is a Lie group and qb is its 
Lie algebra. The Lie bracket in qb maybe expressed in terms of the convolution 
product * in the dual of A (similar developments for the case of words will be 
presented in the next section). 

2.6. The Hamiltonian case. Geometric integration. In many 
applications, the system 0 is Hamiltonian, i.e. the dimension D is even and the 
vector field f(x ) is of the form J~ 1 VH(x ), where H is the (real-valued) Hamilto¬ 
nian function and J is the matrix 


(the four blocks are of size D/2 x D/2). Hamiltonian systems are characterized 
by the property that their solution flow </>ft is, for each h, a canonical or symplectic 
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transformation, i.e. (at each x) 

(d x (j)(x)) T J d x (j>(x) = J. 

It is useful, particularly in the context of long-time simulations, to consider one- 
step integrators x„_|_i = iph{x n ) that when applied to Hamiltonian system originate 
a transformation i^h that is likewise canonical. These integrators are called sym- 
plectic [15], [2Q, [22] and their importance was first highlighted by Feng Kang. It 
was proved independently in ED, [40] . [46] that the RK method © is symplectic 
if 

bidij + bjdji = bibj , 1 < i,j < s. 

In the spirit of the material above, it is sometimes useful to check symplecticness 
by looking at the B-series of the method rather than by examining the method 
coefficients. In m it was proved that, for Hamiltonian systems m, the B-series 
([5|l with Sq = 1 (not necessarily associated with an RK integrator) is symplectic if 
and only if, for each pair of nonempty rooted trees u, v, 

&uov T fi-UOU — (8) 

Here u o v denotes the so-called Butcher product of u and v , i.e. the rooted tree 
obtained by grafting the root of v into the root of u (e.g. . o i = i, : o . = v). In m 
the characterization © was used as a stepping stone to establish the nonexistence 
of symplectic multiderivative RK schemes. 

In terms of modified vector fields, symplectic integrators i [ ph may be charac¬ 
terized as methods that when applied to a Hamiltonian problem m result in a 
Hamiltonian fh ( x ). Thus symplectic integrators may be alternatively seen as those 
integrators that, when applied to a Hamiltonian system, generate a map iph that 
formally coincides with the flow of another Hamiltonian system, hopefully close to 
the true Hamiltonian. For nonsymplectic schemes iph will coincide with the flow 
of a vector field perhaps close to / but not within the Hamiltonian class. This 
interpretation is crucial in symplectic integration )45j , [26] . 

The set of all B-series hfh{x) that are Hamiltonian when / is Hamiltonian 
provides a Lie subalgebra g 0 of qb] the associated Lie subgroup of Qb corresponds 
to of all symplectic iph- An element f3 G qb belongs to go if and only if, for 
nonempty u and v, 

Puov + Pvou - 0 - 

Note that this relation for the Lie algebra is linear, as distinct from the corre¬ 
sponding relation © for the group. For /3 € go the Hamiltonian function Hh(x) 
for the Hamiltonian vector field fh ( x ) is given by a formal series somehow similar 
to © but based on so-called (scalar) elementary Hamiltonians [25] rather than on 
(vector-valued) elementary differentials. While there is an elementary differential 
per rooted tree, there are ‘fewer’ elementary Hamiltonians (one per so-called non- 
superfluous free tree [44]); this explains that for symplectic RK schemes the order 
conditions corresponding to different rooted trees are not independent. 
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The study of symplectic integrators for Hamiltonian problems was the first 
step in what was termed in [41 1 geometric integration: the integration of differen¬ 
tial equations by schemes that preserve important geometric features of the system 
being integrated [26] . In this way, the literature has envisaged volume preserving 
integrators to integrate divergence free differential equations, integrators that pre¬ 
serve relevant invariants, etc. Formal series have been a key element in these 
studies, see e.g. na, na, m- 

2.7. Extensions. There are several useful extensions of the Runge-Kutta for¬ 
mat ([2]). In some problems the components of the vector x come to us divided 
into two or more different groups (for instance in mechanical problems we may 
have positions and velocities). In those circumstances we may use different co¬ 
efficients bi , ciy for different groups of components; the result is a partitioned 
RK scheme. An important particular case arises when a second-order system 
( d 2 /dt 2 )y = F(y, ( d/dt)y ) is rewritten as a first-order system for x = ( y , ( d/dt)y)\ 
this is the realm of Runge-Kutta-Nystrom methods. In other instances the vector 
field in (flj may be decomposed as a sum of N parts 

f(x) = h(x) -1-b f N (x ) 

and we may resort to additive RK methods [T| based on evaluations of the indi¬ 
vidual parts f t rather than on evaluations of /. 

The material outlined in previous paragraphs may be adapted to cover all those 
extensions. A unifying technique has been given in [34] . Let us take the additive 
case as an example. When Taylor expanding the map we find elementary 
differentials like d x fi(x)[f 2 {x)] or d xx f 2 (x)[fi(x) 1 fe^x)}. The structure of these 
is captured by coloured rooted trees, i.e. rooted trees where each vertex has been 
marked with one of the symbols (colours) 1, ..., N (for the elementary differential 
dxxf' 2 (x)[fi(x), f 2 {x)], the root of V is coloured as 2 and the two terminal vertices 
as 1 and 2). In the definition ([5]) one has to replace T by the set of all coloured 
rooted trees; after that, all the developments above are easily adapted to the 
additive case [T. 

In addition to the composition law for B-series that has been the key element 
above, a substitution law has been introduced m , m 

The use of formal series based on rooted trees goes beyond integrators based 
on the RK idea of repeated evaluations of the vector field. Thus the paper [3B] 
studies order conditions for splitting and composition methods. 

To conclude this section we mention the possibility of using the rooted tree 
machinery to perform efficiently high-order averaging in periodically or quasiperi- 
odically forced dynamical systems m, m, ce m- The idea is to expand the 
solution as a B-series with oscillatory coeffients and show that those coefficients 
may be interpolated by nonoscillatory functions of t. As a byproduct one may 
obtain in some circumstances formal conserved quantities. These are other ap¬ 
plications to nonnumerical mathematics of the series methodologies developed to 
analyze numerical integrators. 
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3. Word series 


While the material above dealt with series paratemerized by rooted trees or colour¬ 
ed rooted trees, we now turn our attention to series parameterized by words from 
an alphabet. 


3.1. Definition. In many situations (see [35! . m the problem to be inte¬ 
grated is of the form 

^- x =^2Xa(t)f a {x), x(0)=x 0 , (9) 

q6A 

where A is a finite or infinite countable set of indices and, for each a £ A, X a is 
a scalar-valued function and f a a D-vector-valued map. The simplest example is 
furnished by the autonomous system 

= f a (x) + f b {x), (10) 

where A = {a,b} and A a (t) = Xb(t) = 1. A more complicated nonautonomous 
example will appear in the next section. 

The solution of © has the formal expansion m 


x{t) = x 0 + E E &ai ■■■a n (^)/ai-"a n (^o)) (H) 

n=1 ai,...,a n GA 

where the mappings fa 1 ---aS x ) an d the scalar functions a ai .--a n are defined by the 
recursions 

f ai -a n (x) = d x fa 2 ...a n ( X ) /aifc), 71 > 1, (12) 

and 


ot ai (t) 




) dt\ , 


I ^a n {^n) — dt n ^ Tl 1 . 

Jo 


(13) 


For the particular case m, the inner sum in CD comprises 2” terms and each 
of them has a coefficient t n /n\: the expansion CD is the Taylor series for x(t) as 
a function of t written in terms of the pieces f a and fb rather than in terms of /. 
If the flows </>£ and <j> b h of the split systems ( d/dt)x = / 0 (x) and ( d/dt)x = fb{x) 
are available analytically (or may be easily approximated numerically), it is often 
advantageous to consider integrating (flUl) by means of splitting integrators of the 
form 

i’h = 4>d ah ° ( t , c sh °---o ( t ,b di h ° ( l , c lh ( 14 ) 

with Cj, di, 1 < i < s, constants that specify the method. The Lie-Trotter (jr h o ^ 
and Strang <t>h/2 0 < ° K/2 splittings provide the simplest examples. The Taylor 
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expansion of ifh involves the individual pieces f a and f b and it then makes sense 
to write the Taylor expansion of the solution flow in the form m we have just 
considered. Of course splitting integrators that use the flow of the vector fields f a 
and f b are not to be confused with additive RK schemes that just avail themselves 
of the capability of evaluating the fields f a , f b . 

The notation in dill) may be made slightly more compact by considering A as 
an alphabet and the strings a± ■ • ■ a n as words. Then, if W„ represents the set of 
all words with n letters, CD reads 

OO 

x(t) = X 0 + ^2 5Z aw W fvl X o)- 
n =1 KlGWn 


If we furthermore introduce the empty word 0 and set Wo = {0}, ctg = 1, f®(x) = x, 
then the last expansion becomes 

OO 

x (t) = ^2 53 a ™(t) fw(xo) = Y. ajt) fwjxp), ( 15 ) 

n—0 wGWn w£\V 

where W represents the set of all words. This suggests the following definition: If 
S maps W into C (i.e. <5 € C w )0 then its corresponding word series is the formal 
series 

W&(x) = y S w f w (x). ( 16 ) 

The scalars S w and the functions f w will be called the coefficients of the series 
and the word-basis functions respectively. Thus, for each fixed t, m3 is the word 
series with coefficients a w (t). As we shall see below, in the particular case (flTil) . 
the mapping that represents the splitting method (flTl) corresponds, for each 
fixed h, to a word series. 

The definition of word series is clearly patterned after that of B-series. Note 
however that in 03 the coefficients S w play the role that in ([3 is played by h^S u . 
As we have pointed out already, the expansion 03 corresponds to a family of word 
series: one for each fixed value of t. The reason for this small lack of parallelism 
between the definition of B- and word series is that while © or 0D depend on 
h through powers h J which may be made to feature in the definition, the time 
dependence of (ll5l) and related expansions changes with the functions A a (t). 

Each word-basis function f w , w ^ 0, is build up from partial derivatives of the 
f a , a € A, e.g., if a, b,c G A, 

fba (x) = d x f a {x) f b (x), 

fcba{x) = d x fba(x)f c (x) = d xx f a (x)[f b {x), f c {x)] + d x f a (x) d x f b (x) fc(x). 

Clearly, d x f a (x) f b {x), d xx f a (x)[f b (x),f c (x)], d x f a (x) d x f b {x) f c {x) are elementary 
differentials based on rooted trees coloured by the letters of A. Thus by expanding 

4 Whilc in the case of B-series we only considered real-valued families of coefficients S, it will 
be convenient later to consider word series with complex coefficients. 
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each word-basis function in terms of elementary differentials it would be possible in 
principle to avoid the introduction of words and work with coloured rooted trees. 
However such a move would not be always be advisable because word series are 
more compact than B-series (there are ‘fewer’ word basis functions than elementary 
differentials) and, additionally, the composition rule for word series is simpler than 
its counterpart for B-series. 

3.2. The word series group Q . Given 6,5' £ C w , we associate with them 
its convolution product S * S' £ C w defined by 

n—1 

(5-kS)a = Qj ...q„ + (17) 

3 = 1 

(here it is understood that (5 * 5')® = SqS'q). The convolution product is not 
commutative, but it is associative and has a unit (the element 1 £ C w with 
1 0 = 1 and 1 w = 0 for w / 0 ). 

If w £ W m and w' £ W n are words, m,n > 1, its shuffle product wluw' 
is the formal sum of the (to + n)\/{m\n\) words with to + n letters that may be 
obtained by interleaving the letters of w and w' while preserving the order in which 
the letters appear in each word. In addition 0 luw = w lli 0 = w for each w £ W. 
The operation lli is commutative and associative and has the word 0 as a unit. 
We denote by Q the set of those 7 £ C w that satisfy the so-called shuffle relations: 
70 = 1 and, for each w, w’ £ W, 

N N 

Iwlw’ = Ttoj if wluw' = y^Wj. 

3 =1 3 =1 

The set Q with the operation * is a (non-commutative) formal Lie group, which 
plays here the role played by Qb in the preceding section. For each fixed t , the 
family of coefficients defined by (fl3l) and a@(t) = 1 is an element of the group Q , 
[39l Corollary 3.5]. As we shall see below, some numerical integrators for (0 are 
also associated with families of elements of Q parameterized by the stepsize. 

For 7 £ Q 1 the word series W-y{x) is equivariant with respect to arbitrary (not 
necessarily affine) changes of variables x = x(x) [T7] Proposition 3.1]. In fact, 
if word series are rewritten as B-series (with colours from A ), then the family of 
word series W 1 (x) with 7 £ Q exactly corresponds to the family of B-series that 
are equivariant with respect to arbitrary changes of variables. By implication only 
integrators that are equivariant with respect to arbitrary changes of variables are 
candidates to possess a word-series expansion with coefficients in Q. Thus Q may 
be regarded as a (small) subgroup of Qb- 

In analogy with the composition of B-series, for 7 £ Q, W 1 {x ) may be substi¬ 
tuted in an arbitrary word series Wg(x), 6 £ C w , to get a new word series; more 
precisely 

W s (Wfflx)) =W^s{x), (18) 

i.e. the coefficients of the word series resulting from the substitution are given by 
the convolution product 7 * 5. It follows immediately that, in the particular case 
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m, the map ifh defined in q corresponds, for each fixed h, to a word series 
with coefficients in Q. 

The Lie algebra g of the group Q consists of the elements j3 G C w such that 
/3@ = 0 and, for each pair of nonempty words w, w ', 

N N 

Y 0 Wj =0 if wluw' = Wj. 
j =i j =i 

The bracket operation in g is 


and corresponds to the Jacobi bracket (commutator) of the associated word series, 
i.e. for /?, /?' G g: 

(d x Wp>(x))Wp(x) - (d x Wp{x))Wp>{x) = W[ 0t p'](x). 

Since for 0 G g, the word series Wg (re) belongs to the Lie algebra (for the Jacobi 
bracket) generated by the f a , the Dynkin-Speclit-Wever formula [5S] may be used 
to rewrite the word series in terms of iterated commutators of these mappings: 

oo 1 

n= 1 ai,...,a n Gi4 

(For n = 1 the terms in the inner sum are of the form 0aifai{x).) 

In the particular case where each f a is a Hamiltonian vector field, the iter¬ 
ated commutators are well known to be Hamiltonian and therefore so is Wg (x). 
Furthermore the Hamiltonian function of the vector field Wg ( x ) is 

Ug (x) = Y PwH w (x), 

•uiEVV, 


where, for each nonempty word w = a\ • • • a n , 

H w {x) = -{{• • • {{H ai , H a2 }, H a3 } ■ ■ ■ }, H a }(x). (19) 

n 

Here {■, •} is the Poisson bracket defined by {A,B}{x) = VA(x) T J~ 1 B(x). 

We conclude this subsection with an interpretation of the material above in 
terms of Hopf algebras. The product lli may be extended in a bilinear way from 
words to linear combinations of words. After such an extension, the vector space 
C(A) of all such linear combinations is a unital, commutative, associative algebra, 
the shuffle algebra , denoted by sh(A) (see [39], [35]). Deconcatenation defines a 
coproduct and turns sh(A) into a (commutative, connected, graded) Hopf algebra 
[2j. The sets Q and g are then respectively the group of characters and the Lie 
algebra of infinitesimal characters of the Hopf algebra sh(A). 
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4. Extended word series 

Extended word series, introduced in EH, are a generalization of word series to 
cope with perturbed integrable problems and its discretizations. 


4.1. Definition. We now consider systems of the form 


y 


'O' 

9 


(jj 


f(y,9), 


( 20 ) 


where y G R D ~ d , 0 < d < D, w € is a vector of frequencies ojj > 0, j = 1,..., d, 
and 6 comprises d angles, so that f{y,9) is 27r-periodic in each component of 9 
with Fourier expansion 


f(y, 9) = Y ex P(* k ' 9) Mv) 

kgz d 


(/k(y) and /_k(y) are mutually conjugate, so as to have a real problem). Systems 
of this form appear in many applications, perhaps after a change of variables. For 
instance, any system ( d/dt)z = Mz + F(z), where M is a skew-symmetric constant 
matrix, may be brought to the format (12U1) : other examples are discussed in [57] . 
When / = 0 the system is integrable (the angles rotate with uniform angular 
velocity and y remains constant) and accordingly we refer to problems of the form 
(CH) as perturbed integrable problems and to / as the perturbation (some readers 
may prefer to substitute ef for /). 

After introducing the functions 


/k(2/,0) = exp(*k-0)/ k (y), y G 


D D—d 


, 9 G 


we have 


d 

dt 


-f(v,o) = 

To find the solution with initial conditions 


My, 9). 

kgZ d 


( 21 ) 

( 22 ) 


2/(0) = 2/o, 0(0) = 0o, 

we perform the time-dependent change of variables 9 = rj + tui to get 

d_ 

dt 


(23) 


= Y ex p(* k ■/k(2/, v), 

kgZ d 


a particular instance of ©• The formula (fill) yields 


y(t)' 


2 /( 0 )' 

_y(t)_ 


?/(°) 


+ Y Y «k 1 -k n (0/ki-k n (2/(0),r?(0)), 

n— 1 ki ,...,k n 


where the coefficients a are still given by (fl3l) (with A k (t) = exp(ik • ut)) and the 
word basis functions are defined by m and (fl2l) (the Jacobian in (fl2l) is of course 
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taken with respect to the .D-dimensional variable (y, 9)). We conclude that, in the 
original variables, the solution flow of (l22l) . has the formal expansion 


Mvo,o 0 ) = 


y(t) 


yo 

_|_ 

'O' 



Oo. 

1 

tU! 


£ £ Qkl - kn ( t ) /ki-.-kn (VO, e 0 ). (24) 

n—1 ki,...,k n 


Note that the word basis functions are independent of the frequencies ui and the 
coefficients a are independent of f. 

With the notation of the preceding section, we write (l24l) in the following form 
(here and later x = (y, 9)): 


x{t) = 


fba(i) (Zo) ■ 


We introduce the vector space C = C d x C w and, if {v,S) G C, define its corre- 
sponding extended word series to be the formal series 


T ^ ^ d w f w {x^. 
has the expansion 


^ {y,8) (*^) 

Then the solution (f24|) of (f22l 

z(i) = W (tuMt)) (x 0 ), ( tu,a(t )) € C, 

with a(f) £ G C C w as defined before. 


4.2. The extended word series group Q. The symbol Q denotes the 
subset of C comprising the elements (u, 7 ) with u £ C d and 7 G Q. For each t, the 
solution coefficients ( tu>,a(t )) G C found above provide an example of element of 
Q. Some numerical integrators for (l20l) will be shown below to have, for each value 
of the stepsize h, an expansion in extended word series with coefficients in Q. 

We introduce linear opertors H„, £„ as follows. If v is a d- vector, S„ is the 
linear operator in C w that maps each <5 G C w into the element of C w defined by 
(S„< 5)0 = 6$ and 

(5 v S)w = exp(i(ki H-b k„) -v)S w , 

for w = ki... k„. For the linear operator £ v on C w , (£ v fi)qi = 0, and for each word 
w = ki • • • k„, 

( £,vS)w = *(ki H-b k„) • v S w . 

With the help of E u we define a binary operation if (u, 7 ) G G and (v, S) G C, 
then 

(u, 7)M V , 6) = (l 0 v + S 9 U > 7 * (2„<5)) G C. 

By using m, it is a simple exercise to check that G acts by composition on 
extended word series as follows: 


tT (v,6) (^)) {u, r y)'k(v,5) (*^) 5 


7 &G- 
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In fact we have defined the operation -A- so as to ensure this property. The set Q 
is a group for the product ★ and C d and Q are subgroups of The unit of Q is 
the element 1 = ( 0 , 1 ). 

As a set, the Lie algebra g of the group Q consists of the elements (v,S) G C 
with S G g. For (v, 5), ( u , rf) G C d x g, the Jacobi bracket of the vector fields W( Vt s), 
W may be shown to be given by [37] 

£u8-h8*ri—r)*5) ; 

accordingly the bracket of g has the expression 

[(«> J), (u, tj))] = (0, £ v r) -£, u 5 + 6*r]-r]*6). 

The 0 in the right-hand side reflects the fact that C d is an Abelian subgroup of Q. 

Assume that in (l22l) the dimension D is even with D/2 — d = m > 0 and that 
the vector of variables takes the form 


x=(y,9) = (p 1 ,...,p m - 1 q 1 


i. „ 1 „d. nl 

, a , . . . , CL , C7 , 




where p 3 is the momentum canonically conjugate to the co-ordinate q 3 and a 3 is 
the momentum (action) canonically conjugate to the co-ordinate (angle) 9 3 . If 
each fk(x) in (12111 is a Hamiltonian vector field with Hamiltonian function H^(x), 
then the system (l22l) is itself Hamiltonian for the Hamiltonian function 


d 

J2uja 3 + Hk{x). 
j=i kez d 


For each (w, 0) G g, the extended word series W is a Hamiltonian formal 
vector field, with Hamiltonian function 


d 

J^ujja 3 + fivjHyj , 

j —1 idEW, h ;^0 

with H w (x) as in ©. 

The paper m shows how to use the algebraic rules we have just described to 
reduce ( 1201 ) to normal form where all oscillatory components are removed from the 
solution by means of a suitable change of variables. 


5. Analyzing splitting methods for perturbed inte- 
grable problems 


Splitting algorithms are natural candidates to integrate (l22l) . Given real coeffi¬ 
cients, aj and bj, j = 1,... ,r, we study the splitting integrator 




b {P) 




O 4> a A O ■ 


b {P) 

\h 


O (ft l O 


h {U) 


(25) 


5 Consider the group homomorphism from the additive group C. d to the group of automor¬ 
phisms of Q that maps each p. E into Then Q is the (outer) semidirect product of Q and 
the additive group C d with respect to this homomorphism. 
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where <j)[ U ^ and (j)^ denote respectively the i-Hows of the split systems corre¬ 
sponding to the unperturbed dynamics 


y 


'O' 

e 


UJ 


and the perturbation 


d \y 

Jt e 


f(y,o)- 


Since the unperturbed dynamics with frequencies u>j is reproduced exactly by 
(USD, one would naively hope that the accuracy of the integrator would be dictated 
for the size of / uniformly in w. It is well known that such an expectation is 
unjustified, see e.g. [25] . 1 12] , since the oscillatory character of the solution leads 
to a very complex behaviour of the numerical solution. The algebraic machinery 
of extended word series has been used in 122 to provide a very detailed description 
of the performance of the integrators; we shall quote below some of the results in 
that paper. 


5.1. The truncation error. Clearly, the mapping <f)[ u * has an expansion 
in extended word series 

<t>t U \ x ) = {tu, l)Gg; 

furthermore, 

(j)[ P \x) = W( o, r( t))(x), (0,r(t)) G Q, 

where r(t) G Q comprises the Taylor coefficients, i.e. r w (t) = t n /n\ if w G W„. 
It follows that iph also possesses an expansion as an extended word series with 
coefficients in Q. A simple computation using the definition of ★ shows that: 

W^(hauJ,a(h)) (*£)j 

where a = JT cq, and a(h) G Q is specified by a$(h) = 1 and, for n = 1,2,..., 

5ki - k n {h) =h n ^2 —-— exp(i (cjqki -f- 1 - c jn k„) • uh). 

1 < j i < - - - <j... < r 'h A 


Here, 


and, 


Cj = cii H- 1 - cij , 1 < j < r, 


— j 11 J 1 — • • ' — j n , 

= Vjt+l-jn 11 £ < Tl, jl = ■ ■ ■ = je. < J^+l < ■ ■ ■ < Jn- 

By subtracting the extended word series of the true solution and of the integrator, 
we obtain an expansion of the local error that may be used to investigate the order 
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of consistency m- However such an investigation throws light on the behaviour 
of the numerical method only when \h\ is small with respect to the periods of the 
oscillations in the problem. A much better description of the numerical solution 
may be obtained by means of the modified systems that we describe next. 

5.2. Modified systems. For n = 1,2,..., we look for a system 

= P(h)£9, 

where /3 w (h) = 0 for words with more than n letters, such that, for words with < n 
letters, the extended word series expansion of its flow matches that of the integra¬ 
tor. Such a system may be constructed EZ3 provided that there is no numerical 
resonance of order < n, i.e. there is no set ky, ..., k r with (ki + - • - + k r )-u>h = 2nj, 
r < n, j / 0. Furthermore, in the Hamiltonian case, the modified systems will 
also be Hamiltonian. 

In the limit where n increases indefinitely one obtains, if there is no numerical 
resonance of any order, a modified system whose formal h-flow exactly reproduces 
the extended word series expansion of cj>h- As detailed in m such modified systems 
provide a very accurate description of the behaviour of the computed solutions. 
Among other things, it makes it possible to prove 0{h) error bounds for values of h 
away from first order numerical resonances but not necessarily small with respect 
to the periods of the fastest oscillations in the problem. 
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